Based on the mat1_1_1 and mat1_1_2 loci from Leptographium procerum, I did depth assays to see presence/absence of genes.

After searching for the homolog proteins, I found that the regions of CMW154_Contig94:9951-11159 that includes the two proteins LWAP_08037-RA and LWAP_08038-RA that are the best BLAST hits to the MAT1_1_1 gene and MAT1_2_1

Lets see how they look.

mat1 <- readRDS("mat1_depth.Rds")
mat1.pop2 <- readRDS("mat1_others.Rds")

colnames(mat1) <- gsub(colnames(mat1), pattern = "_dupmrk.bam", replacement = "")

mat1 <- merge(mat1, mat1.pop2, by = c("Chrom","Pos"))

library(reshape2)
library(tidyverse)
## ── Attaching packages ─────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.1     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ tidyr   1.1.0     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(viridis)
## Loading required package: viridisLite
mat1.melt <- melt(mat1, id=c("Chrom","Pos"))

We have 3 genes: MAT1-1-1, MAT1-2-1 and the SLA gene.

Lets see how this looks like:

sla <- mat1[mat1$Pos > 14590 & mat1$Pos < 17898, -c(1,2)]
mat1_1_1 <- mat1[mat1$Pos > 11159 & mat1$Pos < 13083, -c(1,2)]
mat1_2_1 <- mat1[mat1$Pos > 9951 & mat1$Pos < 10800, -c(1,2)]

sla.avg <- apply(sla, 2, function (y){ sum(y > 0) })/nrow(sla) * 100 
sla.avg <- tibble("sample"=names(sla.avg),sla.avg)

mat1_1_1.avg <- apply(mat1_1_1, 2, function (y){ sum(y > 0) })/nrow(mat1_1_1) * 100 
mat1_1_1.avg <- tibble("sample"=names(mat1_1_1.avg),mat1_1_1.avg)

mat1_2_1.avg <- apply(mat1_2_1, 2, function (y){ sum(y > 0) })/nrow(mat1_2_1) * 100
mat1_2_1.avg <- tibble("sample"=names(mat1_2_1.avg),mat1_2_1.avg)     

all.cov <- merge(sla.avg,mat1_1_1.avg, by = "sample" ) %>% merge(.,mat1_2_1.avg, by = "sample" )
all.cov.m <- as.matrix(all.cov[,-1])
rownames(all.cov.m) <- all.cov[,1]

vcfR::heatmap.bp(all.cov.m)

So, in that case:

Mating type 1:

mat1 <- all.cov[all.cov$sla.avg > 99 & all.cov$mat1_1_1.avg > 99 & all.cov$mat1_2_1.avg > 99,]
mat1$mat <- "MAT_1"
mat2 <- all.cov[!(all.cov$sla.avg > 99 & all.cov$mat1_1_1.avg > 99 & all.cov$mat1_2_1.avg > 99),]
mat2$mat <- "MAT_2"

mat.ann <- rbind(mat1,mat2)
mat.ann$mat <- factor(mat.ann$mat)

mat.ann.raw <- mat.ann
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(formattable)

mat.ann <- mat.ann %>% mutate(
    SLA = ifelse(sla.avg > 90,
                  cell_spec(sla.avg, color = "green", bold = T),
                  cell_spec(sla.avg, color = "red", italic = T)),
    MAT1_1_1= ifelse(mat1_1_1.avg > 90,
                  cell_spec(mat1_1_1.avg, color = "green", bold = T),
                  cell_spec(mat1_1_1.avg, color = "red", italic = T)),
    MAT1_2_1 = ifelse(mat1_2_1.avg > 90,
                  cell_spec(mat1_2_1.avg, color = "green", bold = T),
                  cell_spec(mat1_2_1.avg, color = "red", italic = T)),
    mat = color_tile("blue", "orange")(mat)
  )
mat.ann [,-c(2,3,4)]  %>% 
  kable(escape = F) %>%
  kable_styling("hover", full_width = F) %>%
  column_spec(5, width = "3cm")
sample mat SLA MAT1_1_1 MAT1_2_1
CABL1 MAT_1 99.969761112791 99.8439937597504 100
CABL2 MAT_1 100 99.9479979199168 100
CABL3 MAT_1 100 99.9479979199168 100
CADC02 MAT_1 100 99.8439937597504 100
CADC03 MAT_1 99.969761112791 99.5839833593344 100
CADC06 MAT_1 99.9395222255821 99.6359854394176 100
MP01 MAT_1 100 99.9479979199168 100
MP03 MAT_1 100 99.8439937597504 100
MP04 MAT_1 100 100 100
MP06 MAT_1 100 100 100
MP08 MAT_1 100 99.9479979199168 100
MP09 MAT_1 100 99.9479979199168 100
MP10 MAT_1 100 99.8439937597504 100
STF1015-01 MAT_1 100 100 100
STF1015-09 MAT_1 100 100 100
STF1015-13 MAT_1 100 100 100
STF1015-14 MAT_1 99.6068944662836 100 100
STF1015-15 MAT_1 100 100 100
STF1032-07 MAT_1 100 100 100
STF1032-12 MAT_1 100 100 100
STF1032-13 MAT_1 100 100 100
STFO706-01 MAT_1 100 100 100
STFO706-04 MAT_1 100 100 100
STFO706-05 MAT_1 100 100 100
STFO7A2-02 MAT_1 100 100 100
STFO7A2-03 MAT_1 100 100 100
STFO7A2-04 MAT_1 100 100 100
STFO7A2-05 MAT_1 100 100 100
STFO7A2-06 MAT_1 99.6976111279105 100 100
STFO7A2-07 MAT_1 100 100 100
STFO7A2-11 MAT_1 100 100 100
STFO9A5-01 MAT_1 100 100 100
STFO9A5-02 MAT_1 100 100 100
STFO9A5-04 MAT_1 100 99.9479979199168 100
STFO9A5-07 MAT_1 100 99.9479979199168 100
STFO9A5-08 MAT_1 100 100 100
STFO9A5-10 MAT_1 100 99.6359854394176 100
WRTF5180-08 MAT_1 100 100 100
WRTF5180-09 MAT_1 100 100 100
WRTF5180-10 MAT_1 100 100 100
WRTF5180-11 MAT_1 100 99.9479979199168 100
WRTF5180-12 MAT_1 100 100 100
WRTF5360-01 MAT_1 100 100 100
WRTF5360-02 MAT_1 100 100 100
WRTF5360-03 MAT_1 100 100 100
WRTF5360-05 MAT_1 100 100 100
WRTF5360-06 MAT_1 100 100 100
WRTF5360-10 MAT_1 100 100 100
WRTF5360-11 MAT_1 100 100 100
WRTF5360-12 MAT_1 100 100 100
WRTF5480-01 MAT_1 99.969761112791 100 100
WRTF5480-03 MAT_1 100 100 100
WRTF5480-04 MAT_1 100 100 100
WRTF5480-05 MAT_1 99.969761112791 100 100
WRTF5480-07 MAT_1 99.969761112791 100 100
WRTF5480-12 MAT_1 100 99.9479979199168 100
WRTF5481B-02 MAT_1 100 100 100
WRTF5481B-07 MAT_1 100 100 100
WRTF5481B-13 MAT_1 100 100 100
WRTF5481B-15 MAT_1 100 100 100
WRTF8800-02 MAT_1 100 100 100
WRTF8800-03 MAT_1 100 100 100
WRTF8800-05 MAT_1 100 100 100
WRTF8800-06 MAT_1 100 100 100
WRTF8800-09 MAT_1 100 100 100
WRTF8800-12 MAT_1 99.969761112791 100 100
WSTF230-02 MAT_1 100 100 100
WSTF230-03 MAT_1 100 100 100
WSTF230-04 MAT_1 100 100 100
WSTF230-05 MAT_1 100 100 100
WSTF230-06 MAT_1 100 100 100
WSTF230-11 MAT_1 100 100 100
WSTF700-01 MAT_1 100 100 99.8820754716981
WSTF700-03 MAT_1 100 100 100
WSTF700-04 MAT_1 100 100 100
WSTF700-06 MAT_1 100 100 100
WSTF700-07 MAT_1 100 100 99.8820754716981
WSTF700-08 MAT_1 100 99.8959958398336 100
WSTF700-12 MAT_1 100 100 100
WSTF700-13 MAT_1 100 99.9479979199168 100
WSTF7280-01 MAT_1 100 100 100
WSTF7280-02 MAT_1 100 100 100
WSTF7280-03 MAT_1 100 100 100
WSTF7280-04 MAT_1 100 100 100
WSTF7280-05 MAT_1 100 100 100
WSTF7280-06 MAT_1 100 100 100
WSTF7280-07 MAT_1 100 100 100
WSTF7280-08 MAT_1 100 100 100
WSTF7280-09 MAT_1 100 100 100
WSTF7280-10 MAT_1 100 99.9479979199168 100
WSTF7280-11 MAT_1 100 100 100
WSTF7280-12 MAT_1 100 100 100
WSTF7280-13 MAT_1 100 100 100
WSTF730-01 MAT_1 100 100 100
WSTF730-02 MAT_1 100 100 100
WSTF730-03 MAT_1 100 100 100
WSTF730-04 MAT_1 100 100 100
WSTF730-05 MAT_1 100 99.8439937597504 100
WSTF730-06 MAT_1 100 99.7919916796672 100
WSTF730-07 MAT_1 100 100 100
WSTF730-08 MAT_1 100 100 100
WSTF730-09 MAT_1 100 100 100
WSTF730-10 MAT_1 100 100 100
WSTF730-11 MAT_1 100 100 100
WSTF770-01 MAT_1 100 100 100
WSTF770-02 MAT_1 100 99.8959958398336 100
WSTF770-03 MAT_1 100 100 100
WSTF770-05 MAT_1 100 99.9479979199168 100
WSTF770-08 MAT_1 100 100 100
WSTF770-09 MAT_1 100 100 100
WSTF770-10 MAT_1 100 99.9479979199168 100
MP02 MAT_2 100 85.2314092563703 0
MP05 MAT_2 100 85.3354134165367 0
MP11 MAT_2 100 87.1034841393656 0
STF1015-16 MAT_2 100 85.2834113364535 0
STF1015-19 MAT_2 100 85.2834113364535 0
STF1015-2X MAT_2 100 85.2834113364535 0
STF1015-2Y MAT_2 100 85.3354134165367 0
STF1032-01 MAT_2 100 85.2834113364535 24.0566037735849
STF1032-04 MAT_2 100 85.3874154966199 0
STF1032-05 MAT_2 100 87.207488299532 5.54245283018868
STF1032-08 MAT_2 100 87.7795111804472 0
STF1032-09 MAT_2 100 85.2834113364535 0
STF1032-10 MAT_2 100 85.3874154966199 0
STF1032-11 MAT_2 100 85.3354134165367 17.8066037735849
STFO706-06 MAT_2 100 85.2834113364535 0
STFO706-08 MAT_2 100 97.2438897555902 0
STFO706-09 MAT_2 100 85.3874154966199 0
STFO706-11 MAT_2 100 89.0275611024441 0
STFO706-16 MAT_2 100 85.2834113364535 0
STFO706-17 MAT_2 100 85.3874154966199 19.3396226415094
STFO706-18 MAT_2 100 85.3354134165367 0
STFO7A2-01 MAT_2 100 85.3874154966199 0
STFO7A2-08 MAT_2 100 85.3354134165367 0
STFO7A2-09 MAT_2 100 85.2834113364535 35.2594339622642
STFO9A5-06 MAT_2 100 85.2834113364535 0
STFO9A5-09 MAT_2 100 85.2834113364535 20.7547169811321
WRTF5180-01 MAT_2 100 85.2834113364535 0
WRTF5180-02 MAT_2 100 85.2834113364535 0
WRTF5180-03 MAT_2 100 85.2834113364535 0
WRTF5180-04 MAT_2 100 85.3354134165367 0
WRTF5180-06 MAT_2 100 90.3796151846074 9.31603773584906
WRTF5180-07 MAT_2 100 85.3354134165367 13.4433962264151
WRTF5360-07 MAT_2 100 85.3354134165367 0
WRTF5360-08 MAT_2 100 85.2834113364535 0
WRTF5360-09 MAT_2 100 85.2834113364535 0
WRTF5481B-01 MAT_2 100 85.2834113364535 0
WRTF5481B-03 MAT_2 100 85.2834113364535 0
WRTF5481B-05 MAT_2 100 90.275611024441 0
WRTF5481B-09 MAT_2 100 85.3874154966199 0
WRTF5481B-10 MAT_2 100 85.2834113364535 0
WRTF5481B-11 MAT_2 100 85.2834113364535 0
WRTF5481B-12 MAT_2 100 85.1274050962038 0
WRTF5481B-14 MAT_2 100 85.3354134165367 0
WRTF8800-01 MAT_2 100 85.2834113364535 0
WRTF8800-08 MAT_2 100 85.2834113364535 0
WRTF8800-10 MAT_2 100 85.2834113364535 0
WRTF8800-11 MAT_2 98.8811611732688 85.2834113364535 35.6132075471698
WSTF230-09 MAT_2 100 85.2834113364535 27.8301886792453
WSTF230-13 MAT_2 100 93.7077483099324 25.5896226415094
WSTF700-10 MAT_2 100 96.0478419136765 0
WSTF770-04 MAT_2 100 85.2834113364535 0
WSTF770-06 MAT_2 100 85.2314092563703 0
WSTF770-11 MAT_2 100 85.3354134165367 0

Looking at the raw depth of coverage:

MAT1

colnames(mat1.melt) <- c("Chrom","Pos","sample","Depth")
mat1.melt <- merge(mat1.melt, mat.ann.raw[,c(1,5)], by = "sample")

ggplot(mat1.melt[mat1.melt$mat %in% "MAT_1",], aes(x=Pos, y=Depth)) + geom_line() + facet_grid(sample~.) + 
  geom_segment(aes(y=100, yend=100, x=9951, xend=10800, color="MAT1_2_1"), arrow = arrow(length = unit(0.1, "inches")), size=3) + 
  geom_segment(aes(y=100, yend=100,x=11159, xend=13083, color="MAT1_1_1"), arrow = arrow(length = unit(0.1, "inches")), size=3) + 
  geom_segment(aes(y=100, yend=100, x=14590, xend=17898, color="SLA"), arrow = arrow(length = unit(0.1, "inches")), size=3) + theme_classic()

MAT2

ggplot(mat1.melt[mat1.melt$mat %in% "MAT_2",], aes(x=Pos, y=Depth)) + geom_line() + facet_grid(sample ~ .) + 
  geom_segment(aes(y=100, yend=100, x=9951, xend=10800, color="MAT1_2_1"), arrow = arrow(length = unit(0.1, "inches")), size=3) + 
  geom_segment(aes(y=100, yend=100,x=11159, xend=13083, color="MAT1_1_1"), arrow = arrow(length = unit(0.1, "inches")), size=3) + 
  geom_segment(aes(y=100, yend=100, x=14590, xend=17898, color="SLA"), arrow = arrow(length = unit(0.1, "inches")), size=3) + theme_classic()

Proportion of mating types per population

The idea is to test that the mating types are in a 1:1 proportion in the population in order to propose that sexual reproduction is occurring. Lets test that hypothesis:

library(poppr)
## Loading required package: adegenet
## Loading required package: ade4
## Registered S3 method overwritten by 'spdep':
##   method   from
##   plot.mst ape
## 
##    /// adegenet 2.1.3 is loaded ////////////
## 
##    > overview: '?adegenet'
##    > tutorials/doc/questions: 'adegenetWeb()' 
##    > bug reports/feature requests: adegenetIssues()
## Registered S3 method overwritten by 'pegas':
##   method      from
##   print.amova ade4
## This is poppr version 2.8.5. To get started, type package?poppr
## OMP parallel support: available
# Setting up population vectors:
lepto.snp <- readRDS("LeptoPNW.snpclone")

# Removing unused samples
mat.ann.raw <- mat.ann.raw[mat.ann.raw$sample %in% indNames(lepto.snp), ]

Mating proportion per site

pvals.mat <- lapply(levels(lepto.snp@strata$Site), function (x) {
  samples.site <- rownames(lepto.snp@strata[lepto.snp@strata$Site %in% x,])
  mat.sub <- mat.ann.raw[mat.ann.raw$sample %in% samples.site,]
  mat.tab <- table(mat.sub$mat) %>% unname
  max.prop <- max(mat.tab)/min(mat.tab) %>% round(digits = 2)
  expect <- nrow(mat.sub)/2
  chisq.val <- (table(mat.sub$mat)[1] - expect)^2/expect
  chisq.val <- (5 - expect)^2/5
  qchisq(chisq.val, df = 1)
  chisq.val <- (table(mat.sub$mat)[1] - expect)^2/expect + (table(mat.sub$mat)[2] - expect)^2/expect
  chisq.res <- chisq.test(table(mat.sub$mat), p = c(0.5,0.5))
  data.frame("Site"=x, "MAT1"=table(mat.sub$mat)[1] %>% unname, "MAT2"=table(mat.sub$mat)[2] %>% unname, "Proportions"=paste0(round(max.prop,digits = 2),":",1), "Chisq"=chisq.res$p.value)
}) %>% bind_rows()
## Warning in qchisq(chisq.val, df = 1): NaNs produced

## Warning in qchisq(chisq.val, df = 1): NaNs produced

## Warning in qchisq(chisq.val, df = 1): NaNs produced
pvals.mat %>%
  kable(escape = F) %>%
  kable_styling("hover", full_width = F) %>%
  column_spec(4, width = "3cm")
Site MAT1 MAT2 Proportions Chisq
MP 7 3 2.33:1 0.2059032
STF 24 23 1.04:1 0.8840280
WRTF 27 20 1.35:1 0.3072284
WSTF 44 6 7.33:1 0.0000001

Mating proportion per strand

pvals.mat <- lapply(levels(lepto.snp@strata$Stand), function (x) {
  samples.site <- rownames(lepto.snp@strata[lepto.snp@strata$Stand %in% x,])
  mat.sub <- mat.ann.raw[mat.ann.raw$sample %in% samples.site,]
  mat.tab <- table(mat.sub$mat) %>% unname
  max.prop <- max(mat.tab)/min(mat.tab) %>% round(digits = 2)
  expect <- nrow(mat.sub)/2
  chisq.val <- (table(mat.sub$mat)[1] - expect)^2/expect
  chisq.val <- (5 - expect)^2/5
  qchisq(chisq.val, df = 1)
  chisq.val <- (table(mat.sub$mat)[1] - expect)^2/expect + (table(mat.sub$mat)[2] - expect)^2/expect
  chisq.res <- chisq.test(table(mat.sub$mat), p = c(0.5,0.5))
  chisq.res$p.value
  data.frame("Stand"=x, "MAT1"=table(mat.sub$mat)[1] %>% unname, "MAT2"=table(mat.sub$mat)[2] %>% unname, "Proportions"=paste0(round(max.prop,digits = 2),":",1), "Chisq"=chisq.res$p.value)
}) %>% bind_rows()
## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect

## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect

## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect

## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect

## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect

## Warning in chisq.test(table(mat.sub$mat), p = c(0.5, 0.5)): Chi-squared
## approximation may be incorrect
pvals.mat %>%
  kable(escape = F) %>%
  kable_styling("hover", full_width = F) %>%
  column_spec(4, width = "3cm")
Stand MAT1 MAT2 Proportions Chisq
MP_MP 7 3 2.33:1 0.2059032
STF_1015 5 4 1.25:1 0.7388827
STF_1032 3 7 2.33:1 0.2059032
STF_O706 3 7 2.33:1 0.2059032
STF_O7A2 7 3 2.33:1 0.2059032
STF_O9A5 6 2 3:1 0.1572992
WRTF_5180 5 6 1.2:1 0.7630246
WRTF_5360 8 3 2.67:1 0.1316680
WRTF_5480 6 0 Inf:1 0.0143059
WRTF_5481B 3 8 2.67:1 0.1316680
WRTF_8800 5 3 1.67:1 0.4795001
WSTF_230 6 2 3:1 0.1572992
WSTF_700 7 1 7:1 0.0338949
WSTF_7280 13 0 Inf:1 0.0003115
WSTF_730 11 0 Inf:1 0.0009111
WSTF_770 7 3 2.33:1 0.2059032